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ABSTRACT 

The  Digital  Imaging  and  Remote  Sensing  Image  Generation  (DIRSIG)  model  has  been  developed  at  the  Rochester 
Institute  of  Technology  (RIT)  for  over  a  decade.  The  model  is  an  established,  first-principles  based  scene  simulation 
tool  that  has  been  focused  on  passive  multi-  and  hyper-spectral  sensing  from  the  visible  to  long  wave  infrared  (0.4  to  14 
pm).  Leveraging  photon  mapping  techniques  utilized  by  the  computer  graphics  community,  a  first-principles  based 
elastic  Light  Detection  and  Ranging  (LIDAR)  model  was  incorporated  into  the  passive  radiometry  framework  so  that 
the  model  calculates  arbitrary,  time-gated  radiances  reaching  the  sensor  for  both  the  atmospheric  and  topographic 
returns.  The  active  LIDAR  module  handles  a  wide  variety  of  complicated  scene  geometries,  a  diverse  set  of  surface  and 
participating  media  optical  characteristics,  multiple  bounce  and  multiple  scattering  effects,  and  a  flexible  suite  of  sensor 
models.  This  paper  will  present  the  numerical  approaches  employed  to  predict  sensor  reaching  radiances  and 
comparisons  with  analytically  predicted  results.  Representative  data  sets  generated  by  the  DIRSIG  model  for  a 
topographical  LIDAR  will  be  shown.  Additionally,  the  results  from  phenomenological  case  studies  including  standard 
terrain  topography,  forest  canopy  penetration,  and  camouflaged  hard  targets  will  be  presented. 

Keywords:  DIRSIG,  laser  radar  modeling,  ladar,  lidar,  topographic  lidar 


1.  INTRODUCTION 


1.1.  Model  History 

The  initial  development  of  the  Digital  Imaging  and  Remote  Sensing  Image  Generation  (DIRSIG)  model  was  begun  at 
Rochester  Institute  of  Technology  (RIT)  in  the  late  1980’s  as  a  3D  simulation  environment  for  predicting  images  that 
would  be  produced  by  thermal  infrared  systems.  Since  that  time,  the  model  has  been  expanded  to  cover  the  0.35  to  20.0 
micron  region  of  the  spectrum.  The  model  is  designed  to  produce  passive  broad-band,  multi-spectral  and  hyper-spectral 
imagery  through  the  integration  of  a  suite  of  first  principles  based  radiation  propagation  modules.  These  object  oriented 
modules  address  tasks  ranging  from  bi-directional  reflectance  distribution  function  (BRDF)  predictions  of  a  surface,  to 
time  and  material  dependant  surface  temperature  predictions,  to  the  dynamic  viewing  geometry  of  scanning  imaging 
instruments  on  agile  platforms.1  In  addition  to  the  myriad  of  DIRSIG  specific  objects  that  have  been  created,  there  is  a 
suite  of  interface  objects  that  leverage  externally  developed  components  (e.g.  MODTRAN",  FASCODE,  THERM"1) 
that  are  modeling  workhorses  for  the  multi-  and  hyper-spectral  community.  The  software  is  employed  internally  at  RIT 
and  externally  within  the  user  community  as  a  tool  to  aid  in  the  evaluation  of  sensor  designs  and  to  produce  imagery  for 
algorithm  testing  purposes.  Key  components  of  the  model  and  some  aspects  of  the  models  overall  performance  have 
been  gauged  by  several  validation  efforts  over  the  past  decade  of  the  model’s  evolution. IV,V 

1.2.  Historical  Modeling  Approach 

The  modeling  philosophy  that  has  driven  the  DIRSIG  model  development  over  the  years  is  one  that  favors  first- 
principles  radiation  transfer  mechanisms  over  statistical  or  empirical  modeling  approaches.  Most  statistical  or  empirical 
models  have  been  derived  from  specific  data  sets  that  feature  specific  conditions.  The  fit  of  these  models  to  the 


1  Send  correspondence  to  Scott  Brown  at  brown@cis.rit.edu 
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respective  original  data  may  be  exceptional,  but  depending  on  the  underlying  methodology  employed,  the  model  may 
not  be  applicable  at  for  a  different  time  of  day,  for  a  different  season,  for  different  sensor,  etc.  Such  models  may  have  a 
high  degree  of  accuracy  for  specific  cases,  but  very  little  flexibility  for  modeling  others.  These  same  limitations  may  be 
applicable  to  some  classes  of  statistically  based  models.  In  contrast,  the  approach  that  is  employed  by  the  DIRSIG  tool 
is  to  model  as  many  physically  based  interactions  as  possible  by  utilizing  model  inputs  that  predominately  consist  of 
geometric,  optical  and  thermodynamic  inputs.  The  underlying  radiative  transfer  model  then  interacts  with  this  combined 
geometric  and  optical  scene  model  to  predict  the  radiational  flux  into  a  given  direction  for  a  specific  set  of  conditions. 
Although  this  approach  may  not  model  a  specific  data  set  as  accurately  as  an  empirical  or  statistical  that  is  derived  from 
the  imagery,  this  modeling  approach  has  a  higher  degree  of  flexibility  by  allowing  the  user  to  change  the  imaging 
conditions,  scene  conditions,  etc. 

The  general  approach  of  mating  high  fidelity  geometric  models  with  high  fidelity  optical  and  thermodynamic  models 
prevails  throughout  the  DIRSIG  model.  Some  empirically  and  statistically  drive  models  appear  in  the  lower  levels  of  the 
overall  model,  but  the  higher  level  modeling  capabilities  are  derived  from  the  integration  of  a  myriad  of  lower  level 
representations.  The  authors  have  found  this  modeling  architecture  capable  of  reproducing  specifically  sought 
phenomenology  as  well  as  unexpected  collateral  phenomenology  that  might  not  be  realized  using  other  techniques.  The 
images  in  Figure  1  visually  illustrate  some  of  the  spatial  and  spectral  fidelity  resulting  from  the  modeling  approach  used 
by  the  DIRSIG  model  for  an  urban  scene  of  Rochester,  NY. 

1.3.  Active  System  Justification 

Up  until  2002,  the  DIRSIG  model  was  focused  on  simulating  passive  multi-spectral  and  hyper-spectral  sensing  systems 
in  the  visible  through  thermal  infrared  regions.  At  that  time,  the  ability  to  perform  trade  studies  between  passive  and 
active  laser  radar  systems  would  require  the  use  of  separate  modeling  packages  that  may  use  different  geometric,  optical 
and  thermodynamic  descriptions.  In  addition,  the  available  state-of-the-art  LIDAR  scene  simulation  tools  did  not 
support  rigorous  atmospheric  interactions,  participating  mediums,  multiple  bounce/scattering,  thermal  and  reflective 
region  passive  returns,  complex  scene  geometries,  moving  platform  and  scanning  effects,  detailed  material  optical 
descriptions  (BRDF  and  scattering  models)  and  time-gated  returns.  The  benefits  of  an  integrated  active  LIDAR  and 
passive  multi-  and  hyper-spectral  passive  simulation  environment  that  leverage  a  unified  set  of  model  inputs  and 
underlying  radiation  propagation  models  would  be  significant.  A  prototype  of  this  combined  simulation  environment 
was  completed  by  Burton  and  Brown  in  2002. V1,V"  This  paper  discusses  the  fully  developed  model  available  in  the 
current  version  of  the  DIRSIG  scene  simulation  software  with  respect  to  topographic  mapping  applications.  This 
simulation  tool  allows  users  to  simulate  active  and  passive  sensors  using  the  same  scenes  and  scenarios  and  allow  for 
the  evaluation  of  passive  versus  active  approaches  to  specific  problems  and  the  exploration  of  data  fusion  between  these 
two  imaging  paradigms. 


Figure  1:  RGB  images  generated  by  the  DIRSIG  model  of  an  urban  residential  scene  illustrating  the  spatial  and  spectral  fidelity  of 
these  scene  models. 


2.  MODEL  DESCRIPTION 


2.1.  Numerical  Modeling 

The  addition  of  an  active,  laser  radar  capability  to  the  DIRSIG  model  was  accomplished  by  the  addition  of  a  suite  of 
new  objects  to  the  existing  radiometry  framework.  In  general,  the  model  is  designed  to  predict  the  returned  fluxes  from 
the  scene  as  a  function  of  time  with  respect  to  the  shooting  of  the  source  laser.  The  specific  challenges  of  this  imaging 
model  were  driven  by  the  requirement  to  predict  the  received  photon  counts  as  a  function  of  space  and  time.  The  photon 
flux  arriving  at  a  LIDAR  system  often  approaches  discrete  photon  events  due  to  the  low  amount  of  backscattered 
radiation  and  can  prove  difficult  for  many  traditional  Monte-Carlo  ray  tracing  techniques.  Additionally,  the  temporal 
structure  of  these  returns  is  driven  by  the  spatial  structure  of  the  scene  and  the  total  travel  time  of  arriving  photons 
accrued  during  multiple  bounce  and  scattering  events  within  the  scene.  Analytical,  statistical,  and  existing  passive 
radiometry  solvers  were  found  to  be  insufficient  in  many  instances,  particularly  for  low  flux  situations.  Thus,  a  new 
approach  was  sought. 

2.1.1.  Modified  Photon  Mapping 

The  new  approach  that  was  identified  leverages  a  modeling  technique  called  photon  mapping™  Although  the  original 
presentation  focused  on  photon  mapping  as  a  monochromatic  approach,  various  strategies  have  been  introduced  to 
handle  spectral  sources  using  this  methodology.  One  reason  that  this  technique  was  selected  was  because  photon 
mapping  has  been  demonstrated  to  be  applicable  to  traditional  solid  geometiy  reflective  illumination  and  scattering  and 
absorption  by  participating  mediums,  particularly  in  multiple  bounce  and  multiple  scattering  cases. 

The  photon  mapping  approach  is  a  hybrid  of  traditional  forward  and  backward  Monte-Carlo  ray  tracing  techniques. 
In  this  two-pass  method,  source  photons  are  shot  from  a  source  into  the  scene  using  forward  ray  tracing  during  the  first 
pass  and  then  collected  using  a  backward  ray  tracing  during  the  second  pass  (Figure  2).  The  collection  or  rendering 
process  utilizes  the  events  recorded  during  the  first  pass  to  calculate  the  sensor  reaching  radiance.  For  the  purpose  of 
active  laser  radar  applications,  some  modifications  to  the  basic  photon  mapping  treatment  were  made  including  the 
tracking  of  total  travel  time  and  a  literal  photon  counting  process. 


(a)  Pass  1  -  Photon  Tracing  (b)  Recorded  Photon  Map  (c)  Pass  2  -  Photon  Collection 


Figure  2:  Illustration  of  core  photon  mapping  concepts  which  are  (a)  forward  propagation  of  photons,  (b)  the  resulting  photon  map, 
and  (c)  using  photon  map  to  predict  received  photon  counts  at  the  detector. 

During  the  first  pass,  a  modeled  photon  is  cast  into  the  scene  from  the  source  and  performs  a  pseudo-random  walk 
through  the  scene  based  on  the  optical  properties  of  the  surfaces  or  mediums  that  it  encounters.  At  the  location  of  each 
interaction,  information  regarding  the  interaction  event  is  stored  into  a  fast,  3D  data  structure  or  map  ( e,g .  a  kd-treeix,x,xi) 
that  can  be  queried  during  the  second  pass.  A  critical  addition  to  traditional  photon  mapping  paradigm  was  a  field  to 
track  each  photon’s  total  travel  time.  This  time  field  accrues  the  travel  time  for  all  of  the  multiple  bounce/scattering 
events  to  be  recorded  and  is  used  for  time  gating  during  the  second  pass.  The  modeled  photon  is  followed  until  it  is 
absorbed  by  some  element  within  the  scene.  The  photon  casting  process  is  repeated  until  a  specified  number  of 
interaction  events  have  been  recorded  in  the  map.  The  characteristics  of  the  laser  source  are  incorporated  into  the 
spatial,  spectral  and  temporal  distribution  of  the  modeled  photons  fired  into  the  scene. 


From  an  absolute  radiometry  perspective,  each  modeled  photon  represents  a  “photon  bundle”  emitted  from  the 
source.  Much  like  traditional  photon  mapping,  the  variance  of  the  estimated  received  photon  stream  is  intimately  linked 
to  the  number  of  photon  bundles  cast  into  the  scene.  The  number  cast  must  be  sufficient  to  obtain  a  statistically 
significant  number  of  events  throughout  the  field-of- view  of  the  sensor.  Depending  on  the  angular  extent  of  the  beam, 
the  field-of-view  of  the  sensor  and  the  spatial  detail  of  the  scene  approximately  250,000  to  1,000,000  photons  can  be 
utilized  in  the  photon  mapping  process  to  create  reliable  statistics  for  topographic  returns.  For  volume  scattering  returns, 
this  number  may  need  to  be  higher  depending  on  the  absorption  and  scattering  coefficients  and  the  spatial  extent  of  the 
participating  volumes. 

For  the  second  pass,  rays  are  shot  from  the  image  plane  into  the  scene  until  a  scene  element  is  intersected. 
Traditional  photon  mapping  approaches  utilize  the  localized  photon  map  density  to  estimate  the  incident  irradiance  at  a 
surface  or  volume  element.  For  this  application,  the  photon  map  is  used  literally  as  a  photon  counting  mechanism.  The 
model  queries  the  3D  data  structure  that  was  filled  during  the  first  pass  to  find  a  set  of  incident  photons  arriving  at  the 
surface  (or  volume  element  in  the  case  of  a  participating  medium).  Each  of  these  photon  bundles  is  individually 
redirected  towards  the  sensor  using  the  applicable  surface  reflectance  or  volume  scattering  coefficients  for  the  specific 
incident  and  excitant  geometry  (Figure  3).  The  volume  over  which  the  photons  are  collected  is  based  upon  the 
detector’s  field-of-view  and  a  user-defined  spatial  oversampling  factor.  The  received  photon  bundles  are  then 
temporally  distributed  (based  on  the  bundle’s  parametric  temporal  shape)  and  quantized  based  upon  the  user-defined 
listening  window  and  sampling  frequency.  The  result  is  a  time-gated  data  cube  for  the  listening  window  that  contains 
the  received  number  of  photons  at  each  sampling  interval  for  each  detector  element. 


(a)  Surface  (b)  Volume 


Figure  3:  Modified  photon  mapping  photon  count  estimate  and  projections. 

2.2.  Model  Implementation 

The  integration  of  photon  mapping  into  the  DIRSIG  model  entailed  the  implementation  of  several  new  objects.  The 
first  was  the  basic  support  for  photon  mapping.  This  entailed  the  implementation  of  a  3D  data  structure  that  can  be 
quickly  searched  using  spatial  queries  and  the  creation  of  a  photon  object  that  would  be  propagated  and  stored  into  the 
photon  map.  For  this  purpose,  the  traditional  kd-tree  was  used.  The  next  object  was  a  flexible  source  model  that  could 
support  directional  characteristics  and  the  spatial,  spectral  and  temporal  distribution  of  the  source  photons.  In  the  current 
implementation,  the  system  is  modeled  in  a  monochromatic  mode  at  the  peak  wavelength  of  the  source.  The  temporal 
shape  of  the  pulse  is  stored  parametrically  in  each  photon  rather  than  shooting  photons  as  a  function  of  time.  The 
pointing  and  spatial  distribution  of  the  source  is  numerically  modeled  based  on  either  Gaussian  or  top-hat  spatial 
distributions.  The  LIDAR  support  allows  the  user  to  model  co-axial  and  bi-static  systems  by  providing  nearly  arbitrary 
relative  source  to  detector  positioning  and  pointing  geometries. 

2.2.1.  Forward  Propagation 

The  photon  shooting  function  leverages  the  generic  ray  tracing  support  that  already  existed  within  the  DIRSIG  model. 
The  ray  tracer  interacts  with  scene  elements  that  have  material  specific  properties.  Each  material  has  a  set  of  surface 
optical  properties  and  an  optional  set  of  bulk  or  medium  properties.  The  surface  properties  include  a  spectral  reflectance 
and/or  emissivity  property.  The  currently  supported  reflectance  (BRDF)  models  include  importance  based  sampling 
functions  to  support  the  forward  and  backward  Monte-Carlo  ray  tracing.  The  bulk  properties  include  spectral  extinction, 
absorption  and  scattering  coefficient  models.  When  volume  scattering  is  modeled,  a  phase  function  object  can  also  be 
configured  to  describe  the  directional  nature  of  the  scattering.  The  phase  function  objects  also  support  importance 
based  sampling  functions  to  support  the  forward  and  backward  Monte-Carlo  ray  tracing.  Although  these  optical 


descriptions  existed  prior  to  the  addition  of  the  LIDAR  capabilities,  their  interfaces  were  enhanced  during  this  effort  to 
facilitate  an  efficient  implementation  of  the  photon  mapping  subsystem. 

2.2.2.  Backward  Propagation 

The  next  major  component  of  the  model  to  be  implemented  was  the  collection  of  photons  from  the  photon  map  and 
forward  propagation  into  the  imaging  system.  The  modeling  of  the  radiometric  returns  from  a  scene  element  is  handled 
by  a  class  of  objects  in  DIRSIG  referred  to  as  radiometry  solvers.  A  radiometry  solver  encapsulates  an  approach  to 
predicting  the  energy  reflected,  scattered  and  emitted  by  a  surface  or  volume.  For  example,  there  already  existed  a 
specific  solver  for  predicting  the  passive  returns  from  “hard”  surfaces  and  another  for  a  volume  of  gas  ( e.g .  plume). 
These  solvers  utilize  the  material  specific  surface  and  bulk  optical  properties  to  predict  their  results.  One  or  more 
radiometry  solvers  can  be  assigned  to  each  element  in  the  scene.  To  support  the  active  LIDAR  returns,  a  new 
radiometry  solver  was  created  to  compute  the  returns  from  a  scene  surface  or  volume  by  using  the  optical  properties  and 
the  photon  map  to  estimate  the  number  of  incident  photons  at  the  element’s  point  in  space.  Unlike  the  existing  passive 
radiometry  solvers  that  would  place  the  final  result  in  a  time-independent  result  object,  the  LIDAR-specific  radiometry 
solver  places  the  result  into  a  time-gated  result  object.  The  time  gating  is  based  on  a  user-defined  signal  gate  consisting 
of  a  start,  stop  and  delta  time. 

The  final  component  was  the  implementation  of  a  “capture  method”  object  within  the  system  that  directs  how 
outputs  for  the  detectors  on  the  focal  plane  are  computed.  The  LIDAR  specific  capture  method  developed  under  this 
effort  backward  propagates  rays  from  each  detector  element  into  the  scene  where  it  intersects  scene  elements,  runs  the 
appropriate  radiometry  solvers  (passive  and  active),  forward  propagates  the  energy  to  the  focal  plane  and  then  writes  the 
arriving  photons  counts  to  the  output  file. 


2.2.3.  Atmospheric  Returns 

At  this  time  the  atmosphere  is  assumed  to  be  spatial  uniform  as  a  function  of  position  and  altitude,  which  may  not  be 
accurate  for  some  real-world  applications.  The  horizontal  and  vertical  structure  of  the  atmosphere  results  in  different 
absorption  and  scattering  characteristics  as  a  function  of  location.  Ideally,  robust  atmospheric  optical  models  like 
MODTRAN  and  FASCODE  would  drive  both  the  extinction  and  scattering  optical  properties  of  the  atmosphere.  At  this 
time,  the  extinction  coefficients  used  by  the  DIRSIG  model  are  extracted  from  the  existing  MODTRAN  and  FASCODE 
derived  tables.  However,  extraction  of  the  vertically  structured  scattering  coefficients  and  phase  functions  from 
MODTRAN  (FASCODE  does  not  support  scattering)  would  require  custom  modifications  to  the  MODTRAN  code.  For 
the  time  being,  the  user  can  create  their  own  scene  elements  with  user-defined  extinction  and  scattering  properties  to 
replace  the  atmospheric  if  this  level  of  control  is  critical. 

Under  most  conditions,  the  extinction  and  backscatter  coefficients  of  the  atmosphere  are  extremely  small,  which 
means  that  the  probability  of  absorption  and  scattering  events  within  the  atmosphere  is  very  low.  For  example,  the 
scattering  coefficient  for  a  dry  atmosphere  might  be  lxl O’5  [1/m],  which  means  that  you  would  need  to  shoot  105 
photons  into  a  1 -meter  long  box  of  atmosphere  to  get  a  single  scattering  event.  Many  systems  are  attempting  to  resolve 
vertical  resolutions  of  a  fraction  of  a  meter  and  from  an  altitude  of  several  thousand  meters,  which  implies  that  you 
would  need  to  model  10l°  photons  within  each  spatial  detector  element  in  order  to  get  a  single  scattering  event  within 
each  numerical  contribution  element.  To  achieve  robust  statistics,  this  number  would  be  ideally  several  orders  larger. 


To  use  the  numerical  approach  utilized  by  the  photon  mapping  technique,  the  number  of  photons  that  would  need  to 
be  shot  into  the  atmosphere  to  create  a  statistically  accurate  representation  of  the  scattering  events  would  be  many  of 
orders  of  magnitude  larger  than  the  number  of  photons  needed  to  model  the  topographic  returns.  To  avoid  the  problems 
of  predicting  the  atmospheric  returns  numerically,  the  atmospheric  returns  from  the  model  are  currently  modeled 
analytically  using  the  formulation  proposed  by  Measures/11 
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where  P(/lL ,  R)  is  the  total  scattered  laser  power  received  in  watts,  XL  is  the  peak  wavelength  of  the  laser  in  meters;  R  is 
the  range  of  the  contribution  element  in  meters;  PL  is  the  average  power  of  the  laser  pulse  in  watts;  C  is  the  speed  of 
light  in  meters  per  second,  rL  is  the  pulse  length  in  seconds;  A r  /R2  is  the  acceptance  solid  angle  of  the  receiver  optics 
with  a  collecting  area  Ar;  g(XL)  is  the  receiver  optics  transmission  at  the  laser  wavelength;  g(R)  is  the  geometrical 


form  factor  at  range  R;  P(AL,R)  is  the  atmospheric  backscatter  coefficient  at  the  laser  wavelength  and  range  R  in 
inverse  meters  and  cr^ (&L,R)  is  the  atmospheric  total  extinction  coefficient  in  inverse  meters.  For  most  of  the 
topographic  systems  modeled  to  date,  the  atmospherically  scattered  photon  returns  amount  to  only  a  few  photons 
accumulated  over  the  entire  path  length,  which  is  far  below  the  detection  level  of  the  modeled  systems.  However,  if  the 
tools  is  used  to  model  significantly  longer  path  length  or  an  optically  thicker  atmosphere  (containing  fog),  then  these 
numbers  will  grow  to  be  large  enough  for  consideration  by  the  detection  model. 

2.2.4.  Instrument  and  Platform  Modeling 

Most  operational  laser  radar  instruments  are  flown  on  aircraft  and  utilize  some  method  of  aircraft  relative  scanning  to 
increase  the  spatial  coverage  of  the  system.  The  changes  in  viewing  geometry  during  the  scanning  process  and  the 
location,  orientation  and  stability  of  the  instrument  platform  can  affect  the  final  data  products.  For  example,  the  ability 
to  resolve  a  specific  object  in  a  topographic  data  product  derived  from  a  dataset  might  be  dependant  on  the  angle  from 
which  it  is  illuminated.  The  overall  accuracy  of  a  derived  topographic  product  might  depend  on  the  overall  stability  of 
the  platform  and  knowledge  of  the  platform’s  position.  The  DIRSIG  model  has  a  flexible  platform  model  that  allows  the 
platform  to  be  positioned  and  oriented  as  a  function  of  time.  Furthermore,  the  instrument  can  be  pointed  with  respect  to 
the  platform  either  statically  or  dynamically  using  one  of  the  available  instrument  mount  objects.  These  mount  objects 
can  support  temporal  scanning  including  basic  sinusoidal  across-track  scanning  as  a  function  of  a  user-defined  scan  rate. 

2.2.5.  Generalized  Modeling  Process 

With  the  basic  components  of  the  model  now  described,  the  overall  modeling  process  can  now  be  summarized.  A 
modeling  run  consists  of  the  user  specifying  the  scene  to  be  modeled,  the  instrument  and  instrument  mount  description, 
the  source  description,  the  platform  positioning  data,  and  a  set  of  tasks  that  describe  time  windows  over  which  data  is  to 
be  generated.  The  data  generation  process  begins  by  walking  through  each  user-define  task  according  to  a  step  time  that 
is  usually  driven  by  the  source  pulse  rate.  During  each  time  step,  the  platform  and  instrument  mount  positions  and 
orientations  are  computed  for  the  current  time,  the  source  is  fired  which  fills  the  photon  map,  the  focal  plane  is  captured 
which  collects  the  photons  and  propagates  them  to  the  sensor,  and  the  focal  plane  reaching  photon  counts  are  written  to 
the  output  file.  This  cycle  repeats  for  each  time  in  the  task  window  and  for  each  task  in  the  list. 

The  final  product  of  the  DIRSIG  tool  is  a  3D  cube  consisting  of  photon  counts  as  a  function  of  two  horizontal  spatial 
dimensions  and  one  temporal  dimension  (see  Figure  4).  Typically,  the  resulting  data  cube  is  ingested  by  an  external 
focal  plane  model  to  handle  instrument  specific  detection  schemes  ( e.g .  linear  mode  versus  Gieger  mode),  noise 
sources,  etc.  In  these  situations,  the  data  cube  should  be  generated  with  a  significant  amount  over  oversampling  in  the 
spatial  and  temporal  dimensions  to  allow  for  spatial  and  temporal  integration  by  the  external  sensor  model.  Further 
external  processing  of  the  resulting  modeled  raw  instrument  outputs  can  be  used  to  create  final  data  products  (e.g. 
topographic  maps,  gas  detection  maps,  etc.). 

3.  MODEL  DEMONSTRATION  AND  ANALYTICAL  RESULTS 

Although  the  photon  mapping  concept  theory  is  rooted  in  the  underlying  physics  of  the  radiative  transfer  process  as 
outlined  in  Section  2  and  has  produced  visually  appealing  results  for  images  in  the  computer  graphics  community,  the 
question  remains  as  to  whether  the  approach  can  produce  radiometrically  accurate  results,  particularly  for  a  participating 
medium.  The  next  few  subsections  leverage  a  multi-pronged  verification  approach  with  simplified  analytical  scenarios, 
model-to-model  comparisons,  and  sample  phenomenological  case  studies.  One  should  note  that  verification  indicates  a 
proper  physical  basis  such  that  the  results  are  consistent  with  the  underlying  physics. 
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Figure  4:  An  overview  of  an  end-to-end  system  simulation  using  DIRSIG  for  data  generation. 


3.1.  Analytical  Comparisons 


3.1.1.  Simple  Reflective  Surface 


The  governing  equation  for  the  received  signal  captured  by  the  sensor  is  derived  from  the  general  LIDAR  equation  for 
elastic  scattering.  Based  upon  Measures5"1,  the  received  number  of  photons  from  a  range  R  due  to  elastic  backscattered 
radiation  can  be  written  as: 


Ardelec.ed(^)= 
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where  ^,R)  is  a  system  function  determined  by  the  geometric  considerations  of  the  receiver  optics,  the  quantum 
efficiency  of  the  detection  system  at  each  wavelength,  and  the  overlap  between  the  transmitted  laser  beam  and  the  field 
of  view  of  the  receiver;  Ar/R2  is  the  acceptance  solid  angle  of  the  receiver  optics  with  a  collecting  areav4r;  NL{X)  is  the 
average  number  of  photons  in  the  transmitted  pulse  at  wavelength  X;  crex/(A,r)  is  the  extinction  coefficient  of  the 
participating  medium  at  the  wavelength  X  for  the  range  r;  ps(A)  is  the  backscattering  efficiency  at  of  the  target  surface. 
^A,R)  is  usually  assumed  to  be  separable  into  a  wavelength  dependent  weighting;  £(X),  and  a  geometrical  form  factor, 
is  dependent  primarily  upon  the  receiver  design  characteristics  and  for  the  purposes  of  this  subsection  shall  be 
assumed  to  be  unity  whereas  typical  values  range  from  0.5  to  0.8.  The  geometrical  form  factor  is  used  to  adjust  the 
equation  to  account  for  a  variety  of  factors  and  is  often  fairly  difficult  to  evaluate  for  a  real  LIDAR  system.  If  one 
assumes  that  the  laser  power  distribution  is  Gaussian  in  the  target  plane  and  that  the  limiting  aperture  is  the  detector  size 
instead  of  the  telescope  objective  lens  (or  mirror),  then  the  LIDAR  system  function  can  be  expressed  as 
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where/is  the  effective  focal  length  of  the  system,  r0  is  the  radius  of  the  receive  telescope  aperture,  R  is  the  range,  W(R) 
is  the  transmit  beam  waist  radius  at  range  R,  AD  is  the  area  of  the  detector  in  the  focal  plane,  and  A(r,  y/,  r0 )  is  a  circle  of 
radius  r0  that  is  centered  at  (r,  y/)  in  the  target  plane.  When  Equations  2  and  3  are  combined,  the  topographic  LIDAR 
equation  can  be  rewritten  as 
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Using  Equation  4  and  the  system  parameters  described  in  Table  1,  one  can  analytically  evaluate  the  number  of 
photons  detected  from  a  surface.  To  minimize  atmospheric  effects,  a  dry  atmosphere  will  be  assumed  with  the 
extinction  coefficient  of  lxlO'5  [1/m].  The  planar  surface  is  Lambertian  with  a  known  reflectance  and  is  1200  m 
downrange  from  the  co-located  transmitter  and  receiver. 

Table  1.  LIDAR  system  configurations  for  DIRSIG  verification  runs 


Parameter  Configuration  A 


532  nm 


6 


5  KHz 


1  ns 


_ 100  mm _ 

400  mm 


_ 5  mm _ 

8.0  mrad 


Detector  Size  |  [Varied  from  100  pm  to  50  mm’ 

**50  mm  detector  represents  collection  of  entire  beam 


Pulse  Ener 


Pulse  Repetition  Frequency  (PRF 


Pulse  Length  (td) 


Aperture  Radius  (r0)  _ 


Figure  5(a)  displays  both  the  analytical  prediction  dictated  by  Equation  4  and  DIRSIG’s  integrated  output  for  a 
single  pixel  versus  the  surface  reflectance.  As  shown  in  Figure  5(a),  the  analytical  and  numerical  results  very  well 
correlated  and  exhibit  appropriate  behaviors.  Figure  5(b)  is  a  plot  of  the  impact  of  the  detector  size  on  the  number  of 
photons  detected  from  the  Lambertian  surface.  Once  the  detectors  are  collecting  the  majority  of  the  transmit  beam 
(detector  pitch  >  3  mm),  then  the  number  of  photons  detected  varies  linearly  with  the  surface  reflectance  However, 
when  the  detector  sizes  are  such  that  the  receive  fields-of-view  are  small  compared  to  the  transmit  beam  footprint,  the 
geometrical  form  factor  has  a  greater  impact  on  overall  number  of  detect  photons  and  an  incremental  increase  in 
detector  size  results  in  a  dramatic  improvement. 
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Figure  5:  (a)  DIRSIG  integrated  number  of  photons  detected  from  a  Lambertian  surface  of  varying  reflectance  at  1200  m  downrange 
for  a  single  pulse,  (b)  Detector  size  impact  on  the  integrated  number  of  photons  detected  from  a  Lambertian  surface  of  varying 
reflectances  at  1200  m  downrange  for  a  single  pulse.  (Note  that  50  mm  detector  size  used  to  collect  entire  beam,  presuming  a  detector 
limiting  aperture). 


3.2.  Phenomenology  Demonstrations 
3.2.1.  Tree  Crown  Simulation 

The  first  phenomenology  demonstration  will  address  laser  radar  returns  from  a  tree  crown.  For  this  example,  a  single 
tree  is  modeled  and  the  temporal  structure  of  the  reflected  photon  counts  will  be  explored.  The  geometry  of  the  tree  was 
created  using  a  software  package  called  OnyxTREE  and  the  surfaces  were  attributed  with  leaf  spectral  reflectance  and 
transmission  measurements  made  at  RIT.  The  OnyxTREE  software  is  capable  of  producing  physically  realistic  tree 
geometries  using  detailed  descriptions  of  tree  growth  patterns  xm 

Figure  6  includes  a  3D  photon  map  density  plot  resulting  from  a  single  pulse  modeled  with  1,000,000  photon 
bundles.  In  the  plot  of  the  photon  map,  the  shape  and  detail  of  the  tree  crown  can  be  seen  as  well  as  the  resulting 
shadow  on  the  ground.  Within  the  tree  shadow  there  are  some  areas  with  increased  photon  counts  that  are  evidence  of 
direct  illumination  via  foliage  “poke  through”  and  indirect  illumination  due  to  multiple  reflections.  The  plot  of  the 
photon  counts  as  a  function  of  time/distance  in  Figure  6  results  from  spatially  integrating  a  region  encompassing  most 
of  the  tree  (refer  to  the  highlighted  box  in  the  photon  count  frames).  This  magnitude  vs.  time  plot  shows  a  steady  decay 
within  the  tree  crown  due  to  absorption  and  reflection,  the  late  arrival  of  photons  from  ranges  between  the  crown  and 
the  ground  due  to  multiple  bounces  within  the  crown  and  the  ground  return  itself.  The  magnitude  of  the  ground  return 
indicates  a  non-trivial  probability  of  photons  reflecting  off  the  ground  beneath  the  tree  itself.  The  magnitude  of  this 
return  would  be  a  function  of  the  tree’s  optical  properties  and  leaf  density  (leaf  area  index).  The  images  in  the  lower 
portion  of  Figure  6  are  temporal  slices  of  the  photons  counts  arriving  at  the  sensor.  These  spatial  count  density  maps 
reveal  horizontal  and  vertical  structure  of  the  tree  itself  that  can  be  utilized  to  reconstruct  the  tree  height  and  shape.  The 
modeling  of  a  tree  canopy  (instead  of  a  single  tree)  is  included  in  the  end-to-end  scene  simulation  at  the  end  of  this 
paper. 
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Figure  6:  A  demonstration  of  a  single,  deciduous  tree.  The  3D  photon  map  (top  left)  shows  the  density  of  photon  reflection  and 
absorption  effects.  The  plot  captures  the  spatially  integrated  photon  counts  as  a  function  of  time/range  for  an  area  encompassing  the 
tree  (refer  to  the  highlighted  box  in  the  lower  frames).  The  photon  count  frames  (bottom  right)  represent  the  spatial  density  of  the 
photons  arriving  at  the  sensor  for  specific  times/ranges. 


3.2.2.  Camouflaged  Vehicle  Simulation 

Another  common  application  area  for  topographic  laser  radar  systems  is  for  camouflage  penetration  problems.  To 
demonstrate  the  use  of  the  DIRSIG  model  for  this  task,  a  small  scene  was  constructed  that  contains  a  HMMWV  (“hum- 
vee”)  under  a  camouflage  net  held  in  place  by  supports  on  a  terrain  (Figure  7,  top  left).  The  camouflage  net  is  modeled 
as  a  continuous  surface  that  has  holes  cut  into  it  using  a  high  spatial  resolution  “hole  mask”  that  introduces  geometric 
transmission  due  to  spatial  variations  in  fill  factor.  The  solid  areas  of  the  net  are  attributed  with  a  set  of  three  fabric 
materials  that  have  different  reflectance  and  transmission  factors.  The  vehicle  and  surrounding  terrain  is  also  fully 
attributed  with  appropriate  surface  optical  properties. 

The  DIRSIG  model  produced  a  highly  oversampled,  time-gated,  photon  count  cube  for  this  scene.  For  visual 
reference,  the  DIRSIG  model  also  produced  the  height  “truth  map”  shown  in  the  top  right  of  Figure  7  that  illustrates  the 
amount  of  camouflage  “poke  through”  that  can  be  expected  and  the  visibility  of  net  supports  and  vehicle  underneath. 
The  series  of  images  at  the  bottom  of  Figure  7  show  the  spatial  photon  densities  as  a  function  of  time.  At  this  high 
temporal  sampling  rate,  small  differences  in  vertical  structures  can  be  easily  resolved  including  variations  in  the  net 
surface,  the  net  support,  the  roof  and  hood  of  the  vehicle  and  the  terrain  beneath.  The  ability  to  model  scenes  with  these 
complex  interactions  and  surface  properties  accentuates  the  benefit  of  physics  based  approach  over  what  analytical  or 
statistical  models  can  provide. 

3.2.3.  End-To-End  Topographic  Mapping  Demonstration 

The  end-to-end  modeling  of  operational  systems  entails  using  the  DIRSIG  model  to  create  time-gated  photon  counts 
that  can  be  processed  by  an  external  sensor  model  to  produce  simulated  raw  instrument  data.  These  simulated 
instrument  data  sets  can  then  be  used  with  conventional  data  processing  tools  to  create  topographic  data  products.  For 
this  demonstration,  a  data  collection  of  the  ALIRT  system  developed  at  MIT  Lincoln  Laboratories  was  simulated  over  a 
small  scene  located  on  the  RIT  campus.  The  DIRSIG  scene  database  used  in  this  end-to-end  simulation  was  originally 
developed  for  passive,  tower  based  collections  for  the  purpose  of  camouflage  and  landmine  algorithm  testing.  The 
entire  scene  is  approximately  300  x  300  meters  and  contains  terrain,  trees,  and  man-made  elements.  The  operational 
ALIRT  system  was  flown  during  the  winter  of  2004  over  the  corresponding  portion  of  the  RIT  campus  to  capture 
reference  data  sets  for  future  validation  efforts.  The  original  scene  database  was  modified  to  reflect  the  scene  as  it 
appeared  during  the  data  collection  by  the  actual  ALIRT  sensor. 
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Figure  7:  A  demonstration  of  a  HMMWV  under  a  camouflage  net  (top  left).  The  height  truth  map  (top  right)  illustrates  features  that 
are  reproduced  in  the  photon  count  time/range  sequence  (bottom). 
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Figure  8:  The  end-to-end  simulation  of  the  ALIRT  sensor  included  flight  data  from  an  actual  over-flight  of  the  scene  and  the  across- 
track  scanning  of  the  system  (top  left).  The  general  scene  in  the  near-infrared  from  a  tower-based  instrument  (top  right).  The  derived 
topographic  products  viewed  from  overhead  (bottom  left)  and  a  slant  view  (bottom  right).  Note,  that  scene  elements  were  removed 
for  the  active  simulation  to  reflect  the  scene  acquired  by  the  actual  ALIRT  instrument  for  future  validation  efforts.  Topographic 
reconstructions  courtesy  of  ITT  Industries,  Space  Systems  Division. 

The  ALIRT  instrument  uses  a  spatial  array  that  is  scanned  in  the  across-track  direction  from  an  airborne  platform. 
The  across-track  scan  speed  is  slower  than  the  laser  pulse  repetition  rate  so  that  the  collected  area  over  laps  significantly 
from  pulse  to  pulse.  The  ALIRT  instrument  and  collection  characteristics  were  modeled  using  laser  and  instrument 
properties  supplied  by  MIT  Lincoln  Laboratories.  The  simulated  data  cubes  produced  by  the  DIRSIG  model  utilized  the 
flight  data  from  the  actual  ALIRT  over-flight  of  the  scene.  The  resulting  time-gated,  photon  count  cubes  were  then 
processed  using  a  sensor  model  developed  by  ITT  Industries,  Space  Systems  Division,  that  features  a  rigorous  treatment 
of  the  ALIRT  focal  plane  and  instrument.  The  resulting  simulated  raw  instrument  data  streams  were  then  post-processed 
using  QT  Viewer  (developed  at  Johns  Hopkins  Applied  Physics  Laboratory)  to  view  the  resulting  topographic  products 
(see  Figure  8).  The  scan  sweeps  of  the  system  are  clearly  visible  in  the  overhead  and  slant  topographic  projections. 
These  artifacts  arise  due  to  the  noise-free  knowledge  of  the  platform  location  and  platform  relative  pointing  of  the 
instrument  during  scanning.  Future  simulations  will  include  a  noise  term  on  these  quantities  to  introduce  the  inherent 
uncertainty  in  these  values. 


4.  SUMMARY  AND  CONCLUSIONS 

The  extension  of  the  DIRSIG  model  from  passive-only  scene  simulations  to  the  generation  of  active  laser  radar 
simulations  was  presented.  The  numerical  modeling  approach  utilized  by  the  tool  allows  for  the  simulation  of  complex 
scenes  and  in-scene  interactions  that  cannot  be  predicted  using  common  analytical  or  statistical  approaches.  This 
capability  was  achieved  by  extending  the  traditional  photon  mapping  approach  to  track  the  travel  time  of  photons.  The 
model  results  were  compared  to  analytical  predictions  for  specific  source  and  topographic  scene  scenarios  to  verify  the 
absolute  radiometry  of  the  model  implementation.  The  model  was  also  exercised  to  reproduce  common  phenomenology 


including  multiple  bounce  and  foliage  penetration.  Finally,  the  model  was  demonstrated  in  the  context  of  an  end-to-end 
topographic  system  model  by  simulating  the  collection  of  a  real  sensor  (including  platform  motion  and  instrument 
scanning)  over  a  scene  with  significant  spatial  complexity. 

A  rigorous  validation  of  this  DIRSIG  LIDAR  model  for  topographic  applications  is  currently  underway  using  real  - 
world  data  collected  by  the  ALIRT  system.  Other  future  efforts  will  focus  on  the  treatment  of  volume  scattering  and 
absorption  within  this  modeling  approach.  The  support  for  these  radiative  transfer  mechanisms  is  already  included  in  the 
model  and  is  currently  under  evaluation  for  use  in  simulating  Differential  Absorption  LIDAR  (DIAL)  systems.  Future 
modifications  may  include  tracking  the  polarization  state,  the  phase  distributions  and  the  spectral  structure  of  the 
photons  represented  by  each  modeled  photon  stored  in  the  photon  map.  At  this  time,  all  interactions  within  the  model 
are  elastic  in  nature,  but  the  incorporation  of  inelastic  collisions  can  be  handled  within  this  treatment  by  imposing  a 
combination  of  wavelength  and  temporal  shifts  during  each  interaction. 
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